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I. INTRODUCTION 

The dynamics and roughening of moving interfaces in disordered media has been a subject 
of great interest in non-equilibrium statistical physics since the 1980s. Relevant examples 



of physically and technological 



invasion in porous media 



ia fl 3, Q 



y important processes include thin film deposition 1], fluid 



and wetting and propagation of contact lines between 



phase boundaries jf], |jj • The understanding of the underlying physics involved in interface 
roughening is crucial to the control and optimization of these processes. Significant progress 
in the theoretical study of interface dynamics has been made and a number of theories 
have been developed 8J], which in some selected cases agree well with the experimental 
findings i9|. Most of the theoretical understanding in this field is based on modeling interface 
roughening with a local stochastic equation of motion for the single-valued height variable 
of the interface. However, there are several cases of interest where such an approach cannot 
be justified e.g. due to conservation laws in the bulk. This is especially true for processes 
such as fluid invasion in porous media, which is often experimentally studied by Hele-Shaw 
cells 0, II, 12 j or imbibition of paper Q]. It has been shown that in such cases 



spatially local theories cannot provide a complete description of the underlying dynamics. 
For describing the diffusive invasion dynamics in such systems, a phase-field model explicitly 
including the local liquid bulk mass conservation law has been developed and applied to 
the dynamics of ID imbibition fronts in paper ^(|. This was achieved by a generalized 
Cahn-Hilliard equation with suitable boundary conditions, which couple the system to the 
reservoir. Numerical results for roughening from the model are in good agreement with 
relevant experiments 

One of the great advantages of the phase-field approach is that it's possible to analytically 
derive equations of motion for the phase boundaries in the so-called sharp interface limit 
Most recently, we have developed a systematic formalism to derive such equations for 
the 2D meniscus and ID contact line dynamics of fluids in capillary rise [19]. The equations 
are derived from the 3D bulk phase-field formulation, using variational approach as applied 
to relevant Rayleigh dissipation and free energy functionals. Through successive projections, 
equations of motion for the 2D meniscus and ID contact line can be derived. The leading 
terms of such equations (for small amplitude, long wavelength fluctuations) can be shown 
to agree with results obtained from the sharp interface equations in the appropriate limits 



In addition to the need for non-local models to account for mass conservation, in Hele- 
Shaw and imbibition type of problems the inherent quenched disorder should be properly 
taken care of. Unlike thermal disorder, which is relatively easy to handle, quenched disorder 
depends on the height of the ID interface h(x, t) as rj(x, h(x, £)). This makes its influence on 
interface roughening highly nontrivial, often leading to anomalous scaling 21, 2^. Currently, 
for such cases good agreement between theory, simulations and experiments has not been 
achieved. Even on the experimental side some results such as the quantitative values of 
the scaling exponents, are not consistent and difficult to interpret. Very recently, Soriano 
et al. ^i3| conducted an experimental study of forced fluid invasion in a specially designed 
Hele-Shaw cell. The quenched disorder pattern in Hele-Shaw cell is realized by creating large 
number of copper islands that randomly occupy the sites of a square grid on a fiberglass 
substrate fixed on the bottom Hele-Shaw cell. Three different disorder patterns were used. 
Two of them are obtained by random selection of the sites of a square lattice. The third 
kind of disorder is formed by parallel tracks, continuous in the interface growth direction 
and randomly distributed along the lateral direction. It was found that for forced flow the 
temporal growth exponent (3 ~ 0.5 which is nearly independent of experimental parameters 
and disorder patterns. However, the spatial roughness exponent x was found to be sensitive 
to experimental parameters and disorder patterns. Anomalous scaling with \ « 1.0, and a 
local roughness exponent xiocai ~ 0.5 was found in disorder pattern with parallel tracks along 
the growth direction. It was also demonstrated that such anomalous scaling is a consequence 
of different local velocities on the tracks and the coupling in the motion between neighboring 
tracks. 

On the theoretical side, for non-local Hele-Shaw and paper imbibition problems there are 
two different app roaches within the phase-field models to include additive quenched disorder. 
Dube et al. |lfii Ir^ put the quenched disorder inside the chemical potential, the gradient 
of which is the driving force for interface motion. On the other hand, Hernandez- Machado 
* 0i . Q _ d for the effect of ffuetuatioo of He.e-Sffaw g a P thicks as a mo b,fft y 
with quenched disorder in the phase-field model, while keeping the chemical potential free 
of noise. These two approaches lead to quantitatively different roughening properties. 

When considering the problem of capillary rise in a typical Hele-Shaw cell set-up between 
two rough walls more microscopically, the location of the surface of such corrugated walls 
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in the Cartesian coordinate system is a spatially fluctuating quantity, which indicates the 
presence of quenched disorder. An experimental realization is given in Q]. To treat this 
problem faithfully, in solving the phase-field equation such a fluctuating wall surface should 
be treated as a physical boundary without phenomenologically adding quenched noise to the 
equation of motion as done previously. Consequently it is evident that a rigorous analytic 
treatment of such a problem is overwhelmingly difficult. However, in this paper we demon- 
strate that with proper mathematical formulation of the problem, it is possible - albeit with 
some approximations - to analytically derive equations of motion for the meniscus and con- 
tact line dynamics. Most importantly, these equations incorporate the wall disorder in a 
natural way. To achieve this, we utilize an explicit curvilinear coordinate transformation of 
the phase-field equation in order to apply projection methods to unravel the relevant physics 
in the limit of small disorder. To some extent, this kind of curvilinear coordinate transforma- 
tion is similar to the boundary-fitted coordinate system frequently used in Computational 

n 

Fluid Dynamics (CFD) |24|. 

The outline of this paper is as follows: In Chapter|n]we will consider the phenomenological 
2D phase field model of capillary rise with quenched disorder in the mobility, similar to that 



in Refs. 



251 ] . We will adapt the systematic projection method of Kawasaki and Ohta 



to obtain a linearized interface equation (LIE) that describes small fluctuations of an 

nn 

interface, whose deterministic part reduces to the previous result |23l. l25| in a special limit. 
To treat the problem rigorously, In Chapter IIHI we will consider the full 3D phase field 
model with corrugated walls as the source of quenched disorder. The transformation to 
curvilinear coordinates, as discussed above, is introduced to obtain linearized, effective bulk 
disorder from the original curvilinear boundary condition. Following this we develop and 



apply a general projection scheme 



19] to obtain the effective equations of motion for small 



fluctuations of the 2D meniscus, and ultimately for the ID contact line between the meniscus 
and the wall. Again, the deterministic parts of these equations reduce to previously known 
limits in special cases. However, we demonstrate that the forms of the quenched noise terms 
derived here are different from the previous works. 
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II. 2D PHASE FIELD MODEL WITH STOCHASTIC MOBILITY 



re local conservation of bulk mass was in- 
The bulk disorder in their model was 



. 3. 



A 2D phase-field model explicitly including t 
troduced to study capillary rise by Dube et al 
included in the effective chemical potential. Recently a similar model, where the disorder 
was considered through a stochastic mobility coefficient, was studied by Hernandez-Machado 
et al. 13 Ql- I n particular, they assumed a one-sided mobility coefficient, which vanishes 
on one side of the interface. From this model they derived an equation of motion for small 
interface fluctuations. In this section, we will use the systematic projection method intro- 
duced by Kawasaki and Ohta [26], to derive the corresponding linearized interface equation 
(LIE) describing small fluctuations in a sharp interface in a similar model. In our model 

n 

we assume the mobility to_be independent of the phase, as in the previous works 16], but 



spatially stochastic, as in 



23 



251 ] . This corresponds to considering the invading fluid and 
the porous medium, but not the receding fluid. This picture is valid when the receding fluid 
has low density and viscosity. In practice this would mean a gas, such as air, being displaced 
by a liquid, such as water or oil. The model allows a systematic projection of the effective 
noise term at the interface. 

The phase field model describes capillary rise at a coarse-grained level with a phase field 
0(x, t) that obtains the value (f) = — 1 in the phase of the displaced fluid, and <fi = +1 in 
the phase of the displacing fluid. The phase field thus describes the effective component 
densities, and thus must be locally conserved. An energy cost for an interface is included to 
obtain the free energy as 

^[0(x, t)\ = ~[V0(x, t)] 2 + V(0(x, t)), (1) 

where V has two minima at = +1 and = — 1. The details of V are not relevant in the 
sharp interface limit, except to define the surface tension, so we can choose the standard 
Ginzburg-Landau form V{4>) = —<f) 2 /2 + 4 /4 — a<p, where one of the phases can be set 
metastable by nonzero coefficient a. The equation of motion for the conserved phase field is 
given by the continuity equation d t <p = — V- j and Fick's law j = — MV/i, where \x = 8J-/S(j), 
and M = M(l + is the mobility that we choose to be a position dependent stochastic 
variable here. The resulting equation of motion for the phase field is then given by 

d t (j)(x, t) = V • M(x) V/i[0] = MV • (1 + £(x))W [V'(<f>) - V 2 0] , (2) 
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where the variable £ is now the dimensionless, quenched noise. The s harp interface limit of 
this model without the noise is well known, and discussed e.g. in ref. [l8|. The geometry of 
the problem is that of a half-plane, where a reservoir of the displacing fluid is located at the 
x-axis. The boundary condition of the chemical potential at the half-plane boundary can 
be connected to the physical effect that is driving the capillary rise. In this paper we will 
consider spontaneous imbibition, where the rise is driven by a chemical potential difference 
in the medium, which favors the displacing fluid h|. This means that the two minima 
of V are at different heights. In our notation the chemical potential difference is 2a, and 
we consider chemically homogenous medium, where a = const.. Spontaneous imbibition 
corresponds to a Dirichlet boundary condition (/i = const. = 0) at the reservoir [16]. Forced 
flow, where flow is caused by an imposed mass flux into the system from the reservoir, can 
be modeled with the Neumann boundary condition (V/i = Fy), where F is the flux |l7']. An 
analysis along the lines presented in this paper can also be conducted for the case of forced 
imbibition. A recent review of phase field modeling of imbibition is given in Refs. HQ 

Using the Green's function G(r; r') for the 2D Laplacian, equation (J2J) can be inverted 
using Gauss's theorem. This leads to the integro-differential form 

l - J dr'y/d^)Gd t <p' = (1 + - J dr'y/d&WGV'Z' ■ V'// 



- / cZrVdet(p')GyV /2 £' + A. 
Jv 



(3) 



Notation here has been shortened by omitting the function arguments, and using unprimed 
and primed functions for functions of unprimed and primed coordinates, respectively. The 
Green's functions always take both primed and unprimed coordinates as argument. Also 



the coordinate invariant form is used, with integration measure given by a/ det(g). The 
boundary term A vanishes in the case of spontaneous imbibition, or Dirichlet boundary 
condition in half-plane geometry. 

Using the standard ID kink solution method for projection to sharp interface 0,0] in 
normal coordinates gives Eq. Q as 

J dud u (po J ds'du' y/det(g')Gd t (f)' = -(1 + C| u=0 )(o"« + / dud u (p a) 
+ J dud u <Po J ds'du'^det(g')G [d u ^'d u ,^' + (1 - 2u' K')d s ,t,'d s ,K'd u ,<j)' Q ] (4) 
+ / dud u (po [ ds'du'^det(g')GV 2 £(Kd u/ <f)' + a), 
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where the normal coordinates (s, u) are distances along and perpendicular to the interface, 
respectively, k is the local curvature of the interface, a = | J du (d u <j)o(u)) 2 is the surface 
tension of the phase field model, and finally 0o is the ID kink solution d^(po(u) = V'(<po). 
In the Ginzburg-Landau form of V this would be given by (f>o(u) = tanh(-u/v / 2), with the 
appropriate choice of dimensionless units. We have assumed a disorder correlation length 
that is larger than the interface width, which leads to the constant surface tension obtained. 

With two further approximations 3| the standard procedure j^l can be followed. Trans- 
forming the equation to Cartesian coordinates is made somewhat more tedious by the ne- 
cessity to transform derivatives w.r.t. s and u, but standard differential geometry methods 
can be applied. After the sharp interface limit, i.e. 0o — > ~ 1 + 26(w), the transformation 
to Cartesian coordinates, and linearization in small fluctuations of the interface h and the 
noise £, which also eliminates cross terms proportional to h£, we get the LIE as 

J dx' G(x,H ;x',H Q ) + d y G(x,y;x', H )\ y=Ho h(x,t)+ 

d y ,G(x, H Q - x', y') \ y i =Ho h(x', f)] d t [flo(t) + h(x', t)j = (5) 



-ad 2 x h(x,t) -a + -±-^E(x,H (t)), 



where the disorder term H is given by 

2(ar,y) = J dx' J dy'£(x',y')d y/ G(x,y;x',y'). (6) 

Note that the linearization has been carried out in full here. This means that the disorder 
term does not include any dependence on the interface fluctuations. This eliminates the 
non-linearity of the quenched noise, which is one of its characteristic properties, but we 
believe it is not crucial in the regime where the linearization is appropriate. In other words, 
our results show non-trivial features that arise in the effective noise at the interface level 
with this type of multiplicative bulk disorder, even in the linear regime of weak disorder. 

The Green's function for the Dirichlet boundary condition in half-plane geometry is ex- 
plicitly given by 

n( , n 1 , (a: - x'f + (y - y'f 

G{x, y; x , y = — In 7 

47r (x — x'Y + (y + y') z 

Using this, the Fourier space representation of the interface equation (0) becomes 
(1 - e - 2m ° {t) ) d t h(k,t) + \k\d t H (t) (1 + e -2|fc|^oW) h f k t \ = 

(8) 

-a B \k\ 3 h(k,t) + \k\d t H (t)E(k,H (t)) + \k\Ma k , 
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where Ob = Ma, is the Fourier transform of the chemical potential difference (a& = 0, if 
k 7^ 0, when a = const.), and the disorder in Fourier space is given by 

m y) = -\ dy'^k, y') (e-' fc '^') + e-^y-y'^ . (9) 

In the case of columnar disorder, which doesn't depend on y, the interface equation simplifies 
to 

fiMLA \k\dtH (t) (1 + e-Wft)) + <r B \k\\ fu \k\Ma k 
d t h(k,t) = \ _ e -2 mo{t) h ( k > t ) + d t H °(t)£(k)+ 1 _ e _ mHo{ty (10) 

It is noteworthy that the limit k — > the interface equation is 2Ho(t)Ho(t) = ao, readily giv- 
ing the correct Washburn law if we associate lim^o hk(t) = H (t), and lim^o «fc = cto- 



34j Our method of analysis can be applied to the case of forced flow by simply changing the 



boundary condition of the phase field model at the reservoir, and applying the corresponding 
Green's function jl7 |. 

The dispersion relation (JSJ) above is the main result in this section. It involves two length 
scales: a crossover length scale £ x = 2n (^) 2 [16], and the distance from the reservoir Hq. 
The deterministic part of the dispersion relation is plotted in Fig. at the two limits of 
these length scales: The limit H 3> £ x brings out the "deep" limit, kH 3> 1, behavior. The 
limit £ x 3> Hq shows the "shallow" limit, kHo <ti 1, behavior. A plot from the intermediate 
regime with H = £ x is also shown. 

The deterministic part of the dispersion relation here is identical to that previously ob- 
tained by Dube et al. 16( for the case of disorder in the chemical potential. In the "deep" 
limit where d t h = — (^jB\k\ 3 + H \k\j h, our result also reproduces that of Hernandez- 
Machado et al. j3| for the one-sided mobility case. Using different methods the same result 
has also been obtained for the Hele-Shaw setup by Paune and Casademunt llL and Ganesan 

n 

and Brenner |30j[. 

Our disorder term in the LIE is similar to those obtained in Refs. 



11 



23, 301 in the sense 



that in all cases the effective noise is linearly proportional to the velocity of the interface 
propagation. However, the quantitative forms of the noise terms are different when using 
different methods. How these differences influence the kinetic roughening of interfaces would 
need to be determined by extensive numerical comparison between the different results, 
which at this point has not been conducted. As a linear \k\ proportionality in the Fourier 
space representation of the effective noise term is linked to the y-derivative of the Green's 
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FIG. 1: The dispersion relation in Eqs. (|28|) and (|47|) in arbitrary units. The dispersion is 

determined by the two length scales £ x (vertical dashed line) and Hq (vertical dotted line). The 
upper figure focuses on the "deep" regime, with Hq 3> £, x , in the middle figure these length scales 
are the same, and the lowest figure focuses on the "shallow" regime, with £ x 3> Hq. 

function of the Laplacian in real space rep resentation, it appears to us that the linear \k\ 
is present in the noise terms of Refs. and but not in Ref. 2^. The linear \k\ 
dependence is in general characteristic of effective interface noise caused by conserved bulk 
disorder However, the \k\ dependence (\k\d t H E(k, H (t))) dimensionally cancels the 
integral over the kernel in the effective noise, Eq.Q. This is explicit in the case of columnar 
disorder in Eq.(JTUJ), but is equally valid with the noise in the non-columnar case. Thus the 
multiplicative bulk disorder in the mobility leads to a different type of effective noise than 
the chemical potential disorder, which is considered in Q|. Dimensionally this can be seen 
from the definition of the model, Eq. (J2j), where the noise term is in front of the gradient of 
the chemical potential. 

The fact that the columnar disorder leads to effective noise, which is local in Fourier 
space, is in accordance with the concisions of experiments of Soriano * a, Q. and with 
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231 ] . This would indicate that the phase 



the numerical results from the one-sided model 
dependence of the mobility is not crucially important when considering the invasion of a 
viscous fluid into a fluid with negligible viscosity, and when the interface is consequently 
always stable. When the direction of the invasion is reversed, as studied recently with the 

n 

one-sided model in [25], the situation is naturally quite different. 



III. 3D PHASE FIELD MODEL WITH FLUCTUATING WALLS 

While the stochastic mobility case of the previous section is heuristically appealing, a more 
faithful treatment of the wall disorder should start from the microscopic roughness of the 
walls. To this end, in this chapter we will study a 3D version of the same phase-field model 
as the 2D model, but where the mobility is constant and the disorder is explicitly included as 
fluctuations of the wall position. Thus the geometry of the model is that of a Hele-Shaw cell: 
the 3D volume between two walls that are planar on average, but fluctuate. We will show here 
that by proper mathematical formulation this model can indeed be analyzed by a generalized 
projection formalism Q|. The basic idea is to perform a mathematical transformation from 
the basic Cartesian to a local curvilinear coordinate system as defined by the wall itself. To 
this end, we consider the one- wall setup as shown in Fig. 2. The one- wall setup neglects the 
meniscus-mediated interaction between the two contact lines at the two walls. The one-wall 
approximation also neglects finite gap spacing, i.e. the distance between the two walls in the 
Hele-Shaw cell, which fluctuates as result of the wall fluctuations. This induces additional 
disorder effects when the wall fluctuations are comparabale to the gap spacing, but it remains 
to be studied if the gap effect can be separated from the contact line interaction, which would 
be represented by two coupled equations of motion for the two contact lines. In the present 
work, we only consider to the one-wall approximation, or the limit of large gap spacing. 
Disorder at the wall surface is taken into account by describing local corrugations in the 
wall position around y = by a (small) function y = 5H(x,z). The explicit coordinate 
transformation to the local, curvilinear wall coordinate system is defined by 

x' — x , y' = y — 5H(x, z) , z' — z , (11) 

which corresponds to y' = when y = 5H(x,z). This means that in the new coordinate 
system the wall is located back at y' — 0. Given the proper Green's function, Q the phase 
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field equation can be inverted in any geometry and coordinate system as 

| = MVV ^/ d ™<r,n)^ = /*■) + *,, (12) 

where As is the corresponding surface term, and dV\ is the volume element for the coordinate 
system. The Green's function appropriate for the above mentioned coordinate system is 
considered in some detail in Appendix El Here we compute the correction to the original 
Cartesian Green's function to linear order in 5H(x,z). The final result we obtain, after 
neglecting some surface contributions that are discussed in more detail in Appendix El is 
what one would expect by simply plugging the above definitions into the Cartesian Green's 
function and linearizing in 5H: 

G 3D (ri; r 2 ) = G 3D (ri, r 2 ) - 5H(r x )d yi G 3D {ri] r 2 ) - SH^d^G^rr, r 2 ). (13) 

Here G 3 d is the Green's function for the Laplacian in 3D Cartesian coordinates as given 
in Eq. (|A4[) . Here we again only consider spontaneous capillary rise, where the boundary 
conditions for the phase field model are zero chemical potential at the reservoir and zero 
flux at the walls. Thus the surface integral term in Eq.()12|) is identically zero. 



A. Meniscus Dynamics 

The projection and linearization of the integral equation follows the standard projection 



operation theory 



26 



d m the previous chapter for the 2D 
model. The generalization for the present case is straightforward. After projection the 
integral equation is expressed in terms of the 2D meniscus variable H(x, y) and has the 
following form: 

J dx'dy'G m (x, y, H(x, y; t); x', y', H{x\ y'; t))^^- = o b k. (14) 

When linearizing the above equation, it must be done simultaneously in the meniscus fluc- 
tuations, i.e. H(x, y; t) ~ H (t) + h(x, y; t), and in the wall fluctuations using the linearized 
Green's function of Eq. (jl3J) . This results in the linearized Green's function evaluated at the 
meniscus: 

G 3D (x, y, H(x, y);x', y', H(x' } y')) ~ G 3D (x, y, H ; x', y\ H ) 
+5H(x, H )d y G 3D (x, y, H ; x', y', H Q ) + 8H(x', H )d y ,G 3D (x, y, H o] x', y', H ) (15) 
+h(x, y; t)d z G 3D (x, y, z; x', y', H )\ z=Ho + h(x', y'\ t)d z <G 3D {x, y, H ; x', y', z')\ z , =Ho . 
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Substituting this into the meniscus equation (JUj) gives 



I A = <y B d 2 H Q , J B + Ic + Id + Ie + h = cr B X7 2 h(x, y; t) 



(16) 



where the left hand side equation is to the zeroth order, and the right hand side is to the 
first order in h(x,y,t) or SH(x, H n ). These terms are defined in the same fashion as those 



in Cartesian coordinate system 
They are given by 

lA(x,y 

Ic(x,y 
Ii>(x,y 

I F (x,y 



19|. The terms /e and I F arise from the fluctuating wall. 



dxi J dy!G 3D (x,y, Hott)-^!^!, H (t))d t H (t), 
dx 1 / dy 1 d z G 3D (x,y,z;x 1 ,y 1 ,H (t))\ z=Ho h(x,y;t)d t H (t) 



(17) 
(18) 



dxi J dyid Zl G 3D {x,y,H (t)]Xi,yi,zi)\ Zl=Ho h(xi,yut)dtHQ(t) , (19) 



dxi / dy^sD^x^y, H (t);x 1 ,y 1 , H (t))d t h(xi,yi;t) 



(20) 



21^ 



dxi J dy 1 d y G 3D (x,y,H (t);x 1 ,y 1 ,H (t))6H(x,H )d t H (t) , 
d Xl J d yi d m G 3D (x, y, H (t); x u y u H {t))6H(x u H )d t H (t) . (22) 

The zeroth order equation would give the Washburn law, if we used the Green's function 
for the geometry between two walls and assumed a constant curvature for the meniscus. 
We will assume an average profile H (t), which can be considered to obey Washburn's law 
even though we have only a single vertical wall in the system. Since H Q is not needed for 
determining the form of the evolution equation for the fluctuating part h of Eq. (|16J) at that 
single wall, the precise time-dependence of H is not crucial for the analysis to be presented 
below. 

A local equation of motion for the meniscus fluctuations can be obtained by Fourier-cosine 
transformation following Q|. The above terms become 



?*/i*3Wb\ = ^H h(k } t); 
^/k^y% y [Ic} = \e~ 2kH °H h(k,t); 



1 • 



F*I*.3Wd\ = ^h(k,t)(l-e- 2kH °); 



Ie = 0; 

Ho(t) 



2k 



(l-e- 2kHo )5H(k x ,H (t)), 



(23) 
(24) 

(25) 
(26) 
(27) 
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where k = y/k'£ + ky. We then have the meniscus equation of motion using the above in the 

Fourier transform of Eq. (jTSj) : 

^ x kd t H (t) (i + e -^o(t)fc) +a B k 3 - x . - x 

t) = (1 _ e - 2go( t )fc) — Kk, t) + kH 6H(k, H ). (28) 

The deterministic part of the above meniscus equation is identical to the deterministic part 
of the LIE derived from the 2D phase-field model (jSJ), apart from the dimensionality. This is 
by construction, since the same method was used for the same theory in different dimensions 
by applying the corresponding Green's functions. 

A similar analysis can also be performed for the case where the disorder at the walls 
consists of chemical impurities (i.e. spatially fluctuating surface tension) instead of spatial 
roughness In this case the deterministic part of the meniscus equation is by construction 
identical to that of the above. However, there is no effective noise at the meniscus level, 
since the effect of the disorder comes in from the contact line that serves as a boundary 
condition for the meniscus. 



B. Contact Line Dynamics 

To proceed to the level of the ID contact line we consider the generalized variational 
approach Q|. Formally, one can write the 3D phase field model in terms of variations of a 
Rayleigh dissipation functional, and a free energy functional. Then, using approximations 
that express higher dimensional entities in term of the relevant lower dimensional ones, we 
obtain a chain of projection equations as 

5R 3D [<P] 8F 3D [<f>] 



8<j)(x,y,z;t) 6<j>(x,y,z;t) 

5R 2D [H] _ 5F 2D [H] 

* 5H(x,y;t) "' SH{x,y;t) 
_ 5R 1D [C] 8F W [C] 



(29) 
(30) 
(31) 



5C(x;t) 6C(x;t)' 
where R^b refers to the Rayleigh dissipation functional, and F^d refers to the free energy 
functional in d dimensional space. Here the relevant 3D, 2D and ID objects are the phase 
field, the meniscus profile and the contact line profile, respectively. The variable C(x;t) 
denotes the fluctuating contact line profile, and H(x,y;t) = H Q (t) + h(x,y;t) for the one- 
wall case. The corresponding expansion for the contact line is C(x,t) = C (t) + c(x,t). For 
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small fluctuations h and c, consistency requires that C (t) = H (t). The projection from 
3D to 2D is made possible by the ID kink approximation in the direction normal to the 
interface, as demonstrated in the preceding Section. The corresponding approximation we 
have used to make the 2D to ID projection possible is the quasi-stationary (QS) approxi- 
mation V 2 h(x,y;t) = ; h(x,0,t) = c(x,t), which corresponds to the minimum of energy 
constrained by the contact line profile. The meniscus can then be expressed in terms of the 
contact line as 

/OO 
dxi 5f(x-xi,?/)c(x 1 ;t); (32) 
-OO 

g(k x ,y) = e-^ & g{x,y) = -^-j-. (33) 

7r x z — y z 

The explicit forms for the Rayleigh dissipation and free energy functionals that reproduce 
the meniscus equation (JT3)l when plugged into Eq. (jSUj) are 

/oo /*oo P 

dx\dx 2 \ dy\dy 2 \ dtiH{x\,yi,ti)] 
-oo JO J 

xG 3D (xx,yi, H(xi, yx, h);x 2 , y 2 , H(x 2 , y 2 , t))H(x 2 , y 2 , ti) (34) 

/oo poo p 

dx x / dyi / dh ^/l + \VH{x liyi M)\ 2 - (35) 
-oo JO J 

The effective ID functionals can be obtained from the above by inserting the quasi-stationary 
approximation h qs into (|34|) and (|35|) . In order to obtain the ID equation of motion to linear 
order in small fluctuations one needs to expand the functionals to second order in both 
c(x,t) and 5H(x,y), and then take the variation with respect to the contact line as shown 
in Eq.flHH). 

Neglecting the zeroth order equation for the reasons mentioned earlier, the general Fourier 
space equation of motion we obtain for the first order fluctuations is 

Fx/k x [h + h + h + h + h] = -cr B \k x \c(k x , t), (36) 



where the LHS is the variation of the Rayleigh dissipation functional, and the RHS is the 
variation of the free energy. The RHS is recognized as the deterministic restoring force 
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acting on the contact line |32|. The shorthand notations stand for 



I 2 (x) = 2C / dx x dx 2 dx 3 / dy x dy 2 g(x - x x ,y x ) 

J -oo JO 

d Zl G 3D (x x ,y x ,z x ;x 2 ,y 2 ,Co)\ Zl=CQ g{x x - x 3 , yi)c(x 3 , t) (37) 

/CO POO 
dx x dx 2 dx 3 / dy x dy 2 g(x - xi, y x ) 
-oo J0 

d Z2 G 3D (x 1 ,y 1 ,C ;x 2 ,y 2 ,z 2 )\ Z2=Co g(x 2 - x 3 ,y 2 )c(x 3 ,t) (38) 

/oo /»oo 
dx\dx 2 dx 3 I dy x dy 2 g(x — x x ,y x ) 
-oo Jo 

G 3 d(xi, Vi, Co; x 2 , y 2 , C )g(x 2 - x 3 , y 2 )c(x 3j t) (39) 

/CO p CO 

dxidx 2 I dy x dy 2 g(x - x x ,y x ) 
-oo Jo 

d yi G 3D (x x , C ; x 2 , y 2 , C Q )5H(x x , C ) (40) 

/CO p CO 

dx x dx 2 I dy x dy 2 g(x - x x ,y x ) 
-oo Jo 

d y2 G 3D (x x , y x , C ; x 2 , y 2 , C )SH(x 2 , C ). (41) 

Not all of these integrals are solvable in closed form, but can be approximated to a good 
degree of accuracy by the following expressions: 

= ^r-j c (^); ( 42 ) 



2C f°° e -2\kx\C s 

Fx/kAh] = -nrA k ^ / ds * r~? — T 

K\k x \ J x s 3 Vs 2 - 1 



H^ c(fc„t)e- 2 ' 28 IMQ>; (43) 



_ rrl 2c(k x ,t) f°° , l _ e -a|*«|C(b- 

•T^/fc.M = —7^ / ^S- 



^/U*>] = 0; (45) 



!|^) (! _ e -2.28| fe | C0 ) (44) 



^r/fcJ/ 6 ] = -TT-fiobH (k x , C ) / ds _ — =- 
7r\k x \ J x s 2 Vs 2 - 1 

" -C 5H(k x ,C )(l-e- 2 - 5 ^ c °). (46) 



ir\k x 

We note that the corrections to the free energy functional in the curvilinear coordinates are 
of third order in 5H and h. This can be seen by coordinate transforming the area element, 



which in Cartesian coordinates is ^Jl + (Vh) 2 rs 1 + |(V/i) 2 . 

f x _ e -2.B6\k x \c \ 
~ 2 ,V1 "^ (l-e- 2 - 28|MC j 



^^g — 2.S6| fca; | O ^) 

Finally, after approximating « ^ and V — _ 2 , c ^ « 1, the equation of motion for 
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the contact line fluctuations becomes 

3-*\k x \C /, -2.28|fc !t |Co , \ i _ IL 13 o 

= s - K {1 _ e . 2 . 28lk J 0) <k x ,t) + -\k x \C SH(k x ,C ). (47) 

Note that all of the approximations above are for dimensionless quantities, with errors de- 
pending on the physical parameter | A:^ | Co - The relative errors in the approximated functional 
forms are under 3%, when compared against numerical integration of the respective inte- 
grals for different values of |fc x Co|. An exception to this are relative errors of Fx/kJJ^ and 
Fx/kx [Iq] when \k x \ — > 0, as both 1^ and 1$ vanish at the limit, causing the relative error 
behave badly. However, at machine precision away from \k x \ =0 then these errors are no 
more than 15%, and more importantly the error of the complete dispersion relation stays 
within the 3% error margin. This is due to the fact that the dispersion remains finite, as 
one can see from Figure 

Apart from simple numerical factors, the contact line equation above has the same 
functional form as the results derived in the previous sections. In particular, dth = 
— (^Ik-xl 3 + H \k x \ j h in the "ceep" limit k~ l <C H (t), which thus agrees with the pre- 
vious works discussed earlier |h , I221, 0] . This form of dispersion relation is always obtained 
by our method for interfaces in Model B. This has to do with the quasi- stationary ap- 
proximation, which essentially assumes that meniscus fluctuations dampen quickly in the 
direction perpendicular to the contact line, in order to obtain temporally local equations. 
How this leads to the coupling of the meniscus and contact line dynamics is discussed in 
more detail in another publication 2o| . 

The effective noise term we obtain from the 3D model shares the property of linear 
dependence on the velocit y o f the propagation with the 2D mobility noise, and with the 



previous analyses 
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locit y 1 tne propagation witn tne ZD mobility noise, and witn tne 
13, Q^. In the case of surface tension impurities at the wall [3] 
the effective contact line noise is proportional to k^, whereas in the fluctuating wall case 
in Eq. (j47|) the \k x \ dependence is linear. This is analogous to the 2D mobility disorder 
in the sense that the effective noise is different from that obtained for conserved disorder. 
The extra factor of \k x \ as compared to the 2D model (Eq.fjHJ)) comes from the fact that 
the disorder in the 3D model comes from the walls, whereas in the 2D model the disorder 
is in the bulk. We note that the more complicated properties of the noise in the form of 
non-locality in Fourier space are lost by our approximations. Note that the noise is still 
non-local in real space, as is apparent from its real space representation J 6 in Eq. (}4*T]) . In 
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addition to Fourier space non-localities, our one-wall approach doesn't explicitly include the 
gap spacing, which provides an additional physical length scale (^j). In spite of this, our 



results are in accordance with those of Ref. 



Ill ] , where the gap fluctuations were considered 



as the only source of disorder in context of Darcy's law. 



IV. DISCUSSION AND CONCLUSIONS 



In this paper, we have studied the effective interface dynamics of the three-phase con- 
tact line in a Hele-Shaw experiment by deriving the meniscus and contact line equations 
of motion from higher dimensional bulk phase-field theories by projection methods. The 
projection methods take into account the non-local dynamics of the system caused by local 
mass conservation, and can be systematically applied starting from a full 3D description. 
We have here considered two particular models, namely an ad hoc model where the disorder 
is in the effective mobility in 2D 25], and a more microscopic model where wall corruga- 
tion in 3D is explicitly treated with a curvilinear coordinate transformation. In both cases, 
we have focused on the limit of small disorder by linearizing in disorder strength and in 
the fluctuations caused by the disorder. By construction this linearization, performed in 
real space, causes the Fourier space representations of the equations of motion to be local. 
The upside of this is that the effective dynamics are written in a concise manner, and the 
physical predictions are easily interpreted and the equations we obtain are readily affable 
to numerical analysis. The obvious downside is that the procedure involves a number of 
approximations, the validity of which is not certain a priori. 

In particular, the quasi-stationary approximation of Eq. (|32|). which ultimately enables 
our contact line analysis, requires a critical assessment. A more rigorous approach would in 
fact consider the contact line as the boundary condition to the meniscus equation of motion 
(128)1 . However, explicitly solving the meniscus profile as a function of the contact line leads 
to an equation that is, among other things, non-local in time. Thus we are forced to simplify 
the model by using the QS approximation, the validity of which we can consider both from 
a physical perspective, or more rigorously by considering the limits of the meniscus equation 
of motion. Physically, the QS approximation comes from the minimum of meniscus energy 
constrained by the boundary condition of the contact line. This is expected to define the 
meniscus profile when the meniscus moves slowly, and thus it's called the quasi stationary 
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approximation. Mathematically the meniscus equation ([28)1 reduces to the diffusion equation 



when C <^ k 1 <C y a/Co, and d t h(x, y, t) « 0. In this limit the meniscus level disorder 
acts as a source term 

a B V 2 h(x,y;t) = h. (48) 

This leads to an additional disorder term in the QS approximation, which then leads to a 
plethora of new first order disorder terms in Rid and -Fid- However, all these new disorder 
terms arising from Rid are proportional to Cq, and thus not relevant in the QS limit. 
Additionally, the two new disorder terms created in Fid due to Ip cancel each other out 
exactly. Thus, we expect our results with the simplified version of the QS approximation to 
hold in this particular limit. 

In addition to the detailed derivations and new projection formalism presented here, our 
purpose has been to quantitatively compare two different approaches to modeling rough wall 
Hele-Shaw experiments, namely that based on the 2D phase field model with a stochastic 
mobility, Eq.(J2J), and the 3D phase field model with a fluctuating geometry. The projection 
method we use for both cases produces the linear response of the meniscus and contact line 
to small fluctuations. For both cases, the k dependence of the meniscus and contact line 
deterministic LIEs is the same. In particular, in the special case of the "deep" limit where 
k~ l <C H (t), the asymptotic forms of our general dispersion relations are in agreement 
with previous works on the Hele-Shaw problem byPaune and Casademunt Ganesan 
and Brenner and Hernandez-Machado et al. [23]. The main advantage of our method 
is the way it incorporates the noise into the projection, and thus allows us to study the 
effective noise caused to the contact line level by bulk or wall disorder. The main result of 
this analysis is that in both cases the effective noise is linearly proportional to the velocity 
of the interface. While this result qualitatively agrees with the other works cited above, 
quantitative differences remain in the form of the noise terms. The relevance of these 
differences to the actual kinetic roughening of the interfaces remains a challenging numerical 
problem. 
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FIG. 2: A schematic presentation of the curvilinear coordinates considered in the appendix. When 
presented in terms of curvilinear coordinates, the rough wall by definition look straight. However, 
the shift has introduced a coordinate fluctuation in space, where a previously rectangular object 
looks curved when presented in terms of the new coordinates. This gives rise to a bulk repre- 
sentation of the fluctuating wall, which can be analyzed more easily than the original boundary 
condition representation. 
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Appendix 

APPENDIX A: CURVILINEAR COORDINATE SYSTEM BY FLUCTUATING 

WALL 

In this Appendix we will consider in more detail the Green's function of the Laplacian in 
the fluctuating coordinate system 

x' — x , y' = y — 5H(x, z) , z = z . (Al) 

which is schematically depicted in Figure El The generalization to the two-wall setup is 
straightforward, but very tedious including two independent disorder functions. In partic- 
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ular, we will consider the correction to the Green's function to linear order in small wall 
fluctuations 5H(x,z). First, the metric tensor of the above coordinate system can be ob- 
tained by transformation of the Cartesian metric tensor as 

(A2) 

The coordinate transformation from Cartesian coordinates is merely a shift, and thus the 
integration measure should not change, meaning that the Jacobian in the coordinate trans- 
formation of an integral should be unity. This is indeed so, since 

detfl^/j/]) = 1. (A3) 

In the case of Cartesian coordinates we can obtain the Green's function, which we denote 
G3D, using the image charge method with the Dirichlet boundary condition at z = and 
the Neumann boundary condition at y — 0: 

G 3 d = G 3D + G^ D , (A4) 



[9i 



E 



dx 1 dx 1 
dx 1 ' dxi' 



1 + (dJH) 2 dJH dJHdJH 

dJH 1 d z 5H 
d z SH8JH d z 5H l + (d z 5H) 2 



^{x- xi) 2 + (y± yi) 2 + {z- z x ) 2 y/(x- xi) 2 + (y ± y x ) 2 + (z + z x f 



(A5) 

We work in the limit of small fluctuations, so we write the Laplacian in the curvilinear 
coordinates as the Cartesian Laplacian plus a correction, V 2 = V 2 + L. Note that we 
unconventionally denote V 2 = d 2 x +d 2 2 + d 2 3 for any coordinates [xx,X2,Xs\. The correction 
L is explicitly shown below: 



d5H(x, z) 
dx 



d 2 



dxdy 



+ 2 



d5H(x, z) 
dz 



8 2 | d 2 6H(x,z) | d 2 6H(x,z) d 
dzdy dx 2 dz 2 dy 



In order to use Eq. (|12j) for the curvilinear coordinates, we need the Green's function, which 
has the property of [V 2 + L'\ G^j){r\ r^) = —8{r' — r^). Since the Laplacian in the curvilin- 
ear coordinates can be expressed as the Cartesian Laplacian plus a correction, we can find 
the inverse of the curvilinear Laplacian, or G3D, to first order in 5H as 



G-ioir',^) = (V 2 )" 1 - {y' 2 )^ 1 L' (V /2 ) 



(A7) 
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The above operator notation can be written in full form as: 

/oo poo poo 

dx 2 / dy 2 / dz 2 G 3T) (r'; r' 2 )L(r 2 )G m (r 2 ; r[). (A8) 
-oo JO JO 

Substituting G 3 d into the above, we can work out the correction to the Green's function. Af- 
ter using a similar argument to neglect surface integrals of type J dx'G(x, H(x, t); x', O)^(x') 
as we did with Eq.l j l ]) . we find 

G, D {r[\ r' 2 ) = G 3D (r[; r' 2 ) - ^(r^G^; r' 2 ) - 5H(r' 2 )d y ,GMr[; r' 2 ). (A9) 

At this point the primes can just be dropped, since d y = d y >. This result is hardly surprising, 
since a simple substitution of y' = y — 5H(x, z) to G^d{x'\i r 2 ) yields identical results to linear 
order. 

The neglected surface integrals include a reservoir term and a wall term. The reservoir 
term can be readily seen to be small when the meniscus is further away from the reservoir 
than the disorder correlation length. Additionally, the reservoir boundary correction is zero 
when considering columnar i.e. ^-independent disorder. The wall term is more problematic, 
since it involves the boundary correction due to fluctuation in the direction of the wall 
normal. We have to observe the meniscus further away from the wall than the disorder 
correlation length in order to neglect this boundary correction. The absence of boundary 
disorder corrections is highly desirable if we are to keep our formalism tractable, so we have 
neglected the boundary corrections to the Green's function. 
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